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Abstract 

Consider a bunch of interacting electrons confined in a quantum dot. The later is suddenly 
coupled to semi-infinite biased leads at an initial instant t = 0. We identify the dominant 
contribution to the ergodic current in the off-resonant transport regime, in which the discrete 
spectrum of the quantum dot is well separated from the absolutely continuous spectrum of the 
leads. Our approach allows for arbitrary strength of the electron-electron interaction while 
the current is expanded in even powers of the (weak) lead-dot hopping constant r. We provide 
explicit calculations for sequential tunneling and cotunneling contributions to the current. In 
the interacting case it turns out that the cotunneling current depends on the initial many- 
body configuration of the sample, while in the non-interacting case it does not, and coincides 
with the first term in the expansion of the Landauer formula w.r.t r. 

1 Introduction 

The dominant role of electron-electron interaction at mesoscopic scale has long been recognized, 
effects like Coulomb blockade, Kondo correlations or charge sensing being currently observed and 
even manipulated in transport experiments. The typical system consists of a few-level quantum 
dot coupled to source and drain probes (leads). 

In the physics comunity, different approaches to the transport problem in interacting systems 
were developed and intensively used for numerical simulations. The choice of the method depends 
on which parameter of the problem allows for a perturbative treatment. If the interaction strength 
U is rather small, one can use the non-equilibrium Green-Keldysh formalism to compute transient 
or steady-state currents by approximating the interaction effects at different levels [14]. In the 
strongly interacting case, two alternative methods are available: the first one is the T-matrix 
approach [5] , and the second one is the generalized Master equation formalism [TB] . Both methods 
rely on perturbative expansions w.r.t to the parameter r which measures the coupling between 
the leads and the sample, while U is typically much larger that r. 

Compared with the richness of the physics literature on these subjects, only few rigorous results 
exist on time-dependent transport in interacting systems, and they only apply to weakly interacting 
systems [TD1 U\- More precisely, one needs two important conditions in order to guarantee the 
existence of a stationary state (NESS): 1. the single-particle Hamiltonian describing the non- 
interacting system has purely absolutely continuous spectrum and 2. the interaction strength is 
sufficiently small. Under these conditions, one can write down exact formulas for the stationary 
current, but the calculations must be performed perturbatively in the interaction. These results 
are valid both for the partitioning [3. and the partition- free [J] transport scenarios; moreover [TJ, 
one can prove that in the partitioning case the stationary current is independent on the initial 
state of the sample. 

In this paper we deal with a strongly interacting regime which drastically differs from the 
one discussed above. In particular, the previously mentioned two conditions are not satisfied. 
More precisely, here we only consider quantum dots whose discrete spectrum is far away from the 
absolutely continuous spectrum of the leads. Otherwise stated, this off-resonant condition says 
that the single-particle Hamiltonian of the fully coupled system has discrete bound states and no 
resonances. The existence of a stationary state in this case is still an open problem and we do not 
address it here. We know though that even in the non-interacting case one must take the ergodic 
limit in order to kill off the bound state induced current oscillations [TJ [UJ . 
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From now on r will denote the hopping constant between the leads and sample. If I a .t(r) is 
the current at time t > in a given lead a, its ergodic (Cesaro) limit is defined as: 

J a>0O (r) := lim i / J Q , t (r)di f = lim r, [ e-*I a ,t(r)dt) , (1.1) 

where the second equality expresses the known fact that if the Cesaro limit exists, then it can also 
be calculated through the Abel limit r\ \ 0. 

The central object of our study will be the quantity 

/•oo 

J a (»7, r) :=r? / e ~"*/ Q , t (r)^, 77 > 0, (1.2) 



and its behavior as a function of r. We will show that for a fixed 77 > and for r/77 sufficiently 
small one can expand the RHS of Eq. (|1.2[) in a convergent series of even powers of r, that is: 

00 

I a ( V ,r)=Y,r 2k C a Mv)- (1-3) 
fe=i 

The main questions are the following: 

1. When does the Cesaro limit / a)0 o exist? (This would imply that l a (0 + ,r) exists and equals 

la , 00 / • 

2. How many coefficients C ai 2k(v) admit the limit rj \, 0? 

3. If I a , 00 exists, does it have an asymptotic expansion around r = 0? If yes, can we write: 

/a,oo ~ $> 2fc C Q , 2fc (0+)? (1-4) 
fc=l 

We have a good understanding of the problem for non-interacting systems. The off-resonant 
condition on the spectrum of the one particle Hamiltonian is crucial; in the resonant case we can 
prove that I a ,oo exists and has an asymptotic expansion, but its leading coefficient is not given by 
C«,2(0+) (compare (|4.4I) with (I2.15|) V The interacting case is open. 

In this paper we will identify and compute the coefficients corresponding to k — 1,2 for 
the interacting case under some off-resonant conditions. These two terms have a clear physical 
meaning. C a ^. describes the sequential tunneling processes (i.e. electrons enter or leave the dot 
one- by-one). In the off-resonant regime considered here this contribution is absent because the 
energy conservation requires some levels of the isolated dot to be within the continuous spectrum 
of the leads. then gives the dominant contribution to the current and contains the so- 

called cotunneling processes, in which electrons tunnel from and to the dot in pairs (cooperative 
tunneling) . To our best knowledge, the cotunneling regime has not been previously discussed in a 
rigorous context. 

Before summarising the content of the paper let us comment on the different transport regimes 
(i.e. resonant vs. cotunneling) and some subtleties related to the existence of NESS and of the 
perturbative expansion w.r.t. r. In the absence of the electron-electron interaction the steady- 
state current is given by the Landauer formula (see (H3}) which has been rigorously proved using 
various methods [TJ [HJ [BJ H51 03 ■ This formula implies an effective resolvent which is not always 
analytic w.r.t r (see the discussion around (14.2[) and (|4.4p ). 

The paper is organized as follows: in Section [5] we introduce the model, the problem, state the 
main result, and give a number of consequences. Section [3] deals with the thermodynamic limit 
while Section 2] contains the proof of the sequential and cotunneling formulas. Section [S] is devoted 
to numerical results obtained via the generalized Master equation method [T3]. We find out that 
in the off-resonant case the current does not settle to a steady-state, but the ergodic limit seem 
to exist. We conclude in Section |U 
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2 Notation, setting and the main result 



We shall adopt the partitioning approach to the transport problem [3]. A finite system S is 
coupled to M > 2 noninteracting one-dimensional semi-infinite leads (i.e. particle reservoirs) 
at some initial instant to. For simplicity, we consider a discrete model in which the sample 
is modelled as a finite lattice r C 1? and the leads are described by one-dimensional discrete 
Laplacians on the half-line with Dirichlet boundary conditions. The one-particle Hilbcrt space is 
thus H := / 2 (Ni) © • • • © / 2 (N M ) © ^ 2 (r) =: Hl © Us- We shall use the geometrical (standard) 
basis in Hl, which is the set {\i a ) ■ i > 0, 1 < a < M} where i a means the i-th site of the lead 
a. Similarly we have the basis {|m)} me r for 1-Ls- We denote by \m a ) the vector corresponding to 
the sample site to which the lead a is attached. 

The leads are suddenly coupled to the sample at t = 0. Then for t > the single-particle 
Hamiltonian reads as: 

h = h s + h L + h Tl (2.1) 



where 



hs= E t mn \m){n\, (2.2) 

M / \ M 

7 = 1 \i>0 i>l J 7=1 

^T = r^(|0 7 )(m 7 | + |™ 7 )(0 7 |). (2.4) 

7=1 

In the above equation hx is the so called tunneling Hamiltonian and r is the coupling strength. 
Here {t mn } m ,ner is any symmetric matrix and th > is the hopping constant of the leads. 

We also introduce the eigenfunctions and eigenvalues of h$ and the generalized eigenfunctions 
of h 7 : 

h s 4>\ = e\(/)\, h 7 (p J E = E<p^, (2.5) 

where E is the energy associated to an electron propagating on leads with momentum q € (0, n) 
(the leads are identical) . The explicit form of (p^ in a given site i > of the lead 7 is taken to be: 



<pU3)= ,i + 1)q \ E = 2t L cos(q)e[-2t L ,2t L ], |^(0)| 2 = V - ^ ' ■ (2.6) 

When the leads are finite and of length A, the lead spectrum is purely discrete and given by {e 9l } 
where q now takes discrete values. A corresponding eigenfunction is denoted by ip q ^ . The notation 

of the corresponding Hamiltonians is changed into h\ > , and we have: 

^7 A V 9t =£<hW (2-7) 

We now formulate the transport problem in the language of second quantization (see |12j for 
the standard procedures and notations) . Let T = Tl © ?s be the Fock space constructed from 
the Hilbert space H. The interaction of strength U between electrons in the sample is given by 
the two-particle operator: 

V = - E v(m-n)a*{\m))a(\m))a*(\n))a(\n)), (2.8) 

where a*(\m)) and a(|n)) are creation and annihilation operators in the sites m,n and v(m — n) 
is a pair potential which by assumption is bounded for m — n. These operators act in the 
antisymmetrized Fock space JFs generated by l 2 (T). Similarly one defines creation and annihilation 
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operators in the leads, and then rewrites h in the second quantization w.r.t the geometrical basis. 
We use capital letters to denote the second quantized versions of the one particle operators: 
H a = dT(h a ), a e {S,L,T}. Then the total Hamiltonian of the coupled and interacting system 
reads as follows: 

H = H s + V + H L + H T =■ H a + H T (2.9) 

The current operator in the lead a is introduced as the time derivative of the electron number 
operator N a — J2i>o a * (ia)a(i a )- Using the anticommutation relations one gets: 

7 P IP 1 PT 

J a = -eN a = ~[H,N a ] = -j[H T ,N a ] = ^(a* (|O a »a(|m a » - a* (\m a ))a(\O a ))) . (2.10) 

^From now on we adopt the convention e = h = 1 . Note that the same form of J a holds for leads 
of finite length. 

The different chemical potentials of the leads are p, = [pi, p2, •■•> Mm], and the inverse tem- 
perature f} > is taken constant. The equilibrium sub-state of the leads is characterized by the 
following density matrix: 

.(A) ._ tjM 

which consists of a Gibbs state on each lead. 

The initial density matrix of the sample ps can be any positive function of Hs + V, with 
trace one. For example, if at t < the mesoscopic sample is empty, then we have to take 
ps — |0, 0, ...)(0, 0, ...| where |0,0, ...) is the vacuum state in ^5 written w.r.t the occupation 
number basis. But equally well, one may also consider that the sample already contains a few 
interacting particles at t < 0. Let us denote by \v) the eigenstates of Hs + V, and by E v its many- 
body energies {{Hs +V)\v) = E v \v)). Without loss of generality, we will take ps to be a pure 
state given by an initial many-body state (MBS) henceforth denoted by vq. Thus ps — |^o)(^o|- 

The main quantity we are interested in is the statistical average of the current operator on 
lead a. To this end we introduce the statistical operator p^ of the system with finite leads. It 
solves the quantum Liouville equation for t > and is given by: 

^Ht) = e-«» W p^ A \ p^:=p[ A) ®ps. (2.12) 

If B is an observable acting in the Fock space J 7 , we denote by B(t) := e ttH{ ' Be _itfll ' its 
Heisenberg evolution. Then the average value of B at time t is defined as: 



■■= n? = i - afaW „ — (2.11) 



(B{t)) lei := lim Trr{p {A) (t)B} = lim Tr^p^ B{t)}, (2.13) 

A— >oo A— >oo 



whenever this limit exists. Then our results are summarised in the following theorem: 
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Theorem 2.1. Let f(x) — \j{eP x + 1) be the Fermi function and f a (E) = f(E — /i a ). Let xl be 
the characteristic function of the interval [— 2ft,, 2ti,]. Then: 

(i) . The transient current I a ,t(T) in the lead a is given by 

I a ,t{r) := (J a (t)) Ie f, t > 0, (2.14) 
and defines an entire function of r . 

(ii) . Let I a (j),T) = n J Q e~ vt I a ^{T)dt as in (jl.2|) . Then one has (see (|1.3j) and (|2.6|) ): 



C Q , soq :- C Q! 2(0+) - — ^ J 1 - " 1/ ° 4t 2 ^ 

V V ^ 

x {[1 - - ^)]|^ | 2 XL (^ - - - £^)|^J a X£ (^ ~ E VQ )} , (2.15) 
w/jere: 

^, J M = (^,« # M^->, # = *,•■ (2.16) 

(iii). Assume that the following two off-resonant conditions are fullfiled: 

a) . If \v), \v') differ by one particle, then E v — E u > ^ [— 2tL, 2ti\; 

b) . If \u), \v') differ by two particles, then E v — E v * ^ [— \t^,\ti\. 
Then we have: 

1 f 2th ( E 2 \ 

C a , cot := C Q , 4 (0 + ) = -271- E / dE 1 - W) {V ^ {E) ~ ' ( 2 ' 17 ) 

where V 1QL is the cotunneling rate: 

V ia (E)= Y, ( 2 - 18 ) 
it? t? tp \ f~t( E )\ l ~ fai E - E w + E vo )] 

jj/jj// [771 j i/q [77L a ) 



■r i /- r > t 1 - fi{E)]f a {E + E v > - E v ) . 
-XL{E + E y , - E Uo ) — _ ^ + ff)^, _ ^„ + £; ) ^o^( TO 7)^'( TO «)^V"(rrt 7 )A t/ , Vo (m a ) f> . 



Remark 2.2. Provided that (|1.4I) holds true, if r is sufficiently small then the ergodic current I ai00 
should be well approximated by r 2 C Q , i 2(0+) + r 4 C Q , i 4(0+). Both terms describe tunneling processes 
from and into the dot. Note that E v > = E Vo is allowed in the above sums, thus xl(E — E v i + E Uo ) 
has to be replaced by 1 in those terms. 

Remark 2.3. In the expression of C aj 2(0+), the factor (1 — f a (E))\A VVo \ 2 is the tunneling prob- 
ability from the dot to the leads of an electron with energy E (the corresponding state in the 
lead must be empty). Similarly, the second term of (7,^2(0+) represents processes in which the 
many-body state of the dot changes by 'absorbing' one electron from the leads. These processes 
are called sequential, as electrons tunnel one by one. It is clear that in the off-resonant regime 
(i.e. E v — E Ua ^ [— 2tL : 2ti]) the sequential tunneling is suppressed and one has to go to the next 
term. Note that in the resonant regime of the non-interacting case, this term cannot be recovered 
by expanding the Landauer formula in powers of r (see (|4.4j) for further details) . The phenomenon 
which happens is well described by the following toy example. Let K be a constant either equal 
to or 1. Define the functions: 

?7r 2 r 4 / 1 \ 

I(n, r; K) := (1 - K)-!—^ + 2 arctan 2 . 

n + r A K+n + T A \K+rj + T*J 
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The resonant case is modeled by the condition K — 0. In that case we have: 

7(0+, r;0) = r 2 arctan(l/r 2 ) =r 2 |+C(r 4 ), i"(r?,r;0) =t 2 +C(t 4 ), 

which shows that C aj2 (0+) = 1 and we cannot recover the 'true' behavior of J(0 + ,r;0) from such 
an expansion. 

The off-resonant case is modeled by K = 1. Then: 

7(0+, r; 1) = r 4 - + 0(r 6 ), 7(r?, r; 1) = -I— arctan f -L-) + 0(r e ). 
4 I + 77 \l + r)J 

In this case we see that C Qi 2(0+) = and C a ^(Q+) = 7r/4, and they provide a good approximation 
for the 'true' value of 7(0+, r; 1). 

Remark 2.4. The contribution 7 Q!CO t := r 4 C Q . 4(0+) is the so-called cotunneling current. Further 
discussion on it will be given in Section 5.2. Here we only stress that in the absence of the bias, 
7 Q ,cot = because in this case the chemical potentials of the leads are equal and hence V ia = 7 J, a7 . 
In the off-resonant non-interacting case, we can prove that it docs not depend on the initial state 
in the sample and it is given by the first term of the Landauer formula (see (14.171) ). 

Remark 2.5. Memory effects and dependence on the initial state. For small samples one 
is able to simplify the formula giving the cotunneling current. A typical example is a two-site 
quantum dot. Let us denote by e\ y2 the eigenvalues of the non-interacting dot. We also have that 
hs4>i — e i0i an d hs<t>2 = e202- The four many-body states are E\ = (empty sample), E2 — e\ 
(the ground state of hs), E 3 = e2 (the excited state of hs) and E4 = e\ + e 2 + U (fully occupied) 
where U denotes the strenght of the Coulomb interaction. Let us consider that the initial state 
of the system is |^o) = 1 10) and E Vo = e\, which means that before the coupling we start with 
exactly one electron in the sample, occupying the lowest level. 

The two spectral conditions imposed by the off-resonant regime have to be checked for any 
given set of parameters. Let us explicitely write down these conditions for a sample having only 
two sites: a). e 1>2 i [-2t £ ,2t £ ], e lt2 + U $ [-2t L ,2t L ]; b). E 4 - E x = a + e 2 + U $ [-At L ,4t L ]. 
These conditions can be satisfied in many situations, for example when both are either very 
negative or very positive such that they are far away from the spectrum of the leads. 

The cotunneling current in Eq. (|2.18[) can be further simplified by calculating the coefficients 
A and A* . In order to do that we have to express the creation and annihilation operators in the 
contact sites a#(\m a )) and a^(|m 7 )) in terms of creation and annihilation operators in a given 
single-particle eigenstates a#(\(f>i)) and a#(|</>2)). This leads to obvious selection rules for the 
MBS. Calculating the cotunneling rates terms one can identify elastic and inelastic contributions 
to the current: 

7 Q ,cot = t 4 C q ,4(0+) = 7 e i + 7 in , (2.19) 

where: 



7 R , = 



E 



2L L 



'21 [_ 



dE 1- 



E-ei 



E-e 2 -U 



(ME) - ME)) 



(2.20) 



and 
7 m 



21 L 



dE 1 



XL(E + e! -e 2 0i(m 7 \^ 2 {m a )\ - 1 — 

e 2 + U — E 

1 Tp I \| 1 / M2 U / x, 2 /a( £; + e 2-ei)(l-/ 7 (7;)) 

Xl(E + e 2 - ei)|0i(m 7 )| \<p 2 {m a )\ 



1 



1 



e 2 - E e 2 + U-E 
1 1 



d-E 



e x + U-E ei-E 



{a <-> 7}. 



(2.21) 



G 



Let us comment on the two contributions to the cotunneling in this case. Obviously / i is given 
by a Landauer formula, even if the interaction strength U appears in one of the denominators. 
The two electrons implied in the pairwise tunneling have the same energy E hence this is elastic 
cotunneling. This contribution can be compared with the one calculated in Ref. [16) via what 
the authors call the 'T-matrix method'. To make the connection to their results one should use 
the cotunneling rate given in Eq. (19) of Ref. [16] and calculate the steady-state current as 

111 111 ■ 

In contrast, I m can no longer be written in a Landauer form and contains inelastic processes, as 
the energies in the Fermi functions do not coincide. Note that I ln vanishes in the non-interacting 
case: this happens because for U — the contributions of various inelastic processes cancel each 
other. Moreover, if |e2 — e\\ > At^ then xl{E + e-i — e\) and xl{E + e\ — e?) will vanish for all 
E e [— 2t^, 2ti], thus again ij n = 0. But otherwise it is nonzero. 

We can repeat this computation choosing the initial condition |fo) = 1 00) (the sample is empty 
before coupling it to the leads). In this case we find: 

2 

ii(m a )0i(m 7 ) 4>2(m a )(j)2(m 7 ) 



2t L / p2 



Ia ' cot - L L dE v * 



(UE) - ME))- 



E — e\ E — 62 

(2.22) 

Otherwise stated, for this initial state of the sample the cotunneling current is given by the non- 
interacting Landauer formula. This means that the cotunneling current in the interacting case 
depends on the initial conditions of the sample. This is not such an unexpected result, as different 
initial many-body configurations of the sample select different relevant cotunneling processes. We 
stress though that this memory effect concerns only the cotunneling current in the off-resonant 
regime. In the resonant case where sequential and cotunneling processes coexist we do not expect 
this to happen. 



3 Proof of (i): thermodynamic limit and the definition of 
the transient 

In mesoscopic quantum transport we have to deal with two aparently contradictory conditions: 
1). the leads must be finite if we want the total density matrix to be trace class, and in that 
case the total Hamiltonian has purely discrete spectrum; 2). the total Hamiltonian must also 
have some continuous spectrum since otherwise the ergodic current would be identically zero. The 
correct way out is to fix the time t, define the expectations at finite leads and afterwards make 
them infinitely long. Only after the thermodynamic limit we can let t to go to infinity. More than 
that, the total density matrix is not the good object to work with, and any formal perturbative 
expansions in r at t — oo before the thermodynamic limit has no clear mathematical meaning. 

In this section unless otherwise stated the leads are assumed to be of finite length A. But 
for the simplicity of writing we omit the label A on the leads' Hamiltonian. In order to get an 
expansion of the current in powers of the tunneling Hamiltonian we define W(t) = e e~ , 
verifying the equation: 

iW(t) = H T (t)W(t), W(G) = 1, H T (t) := e itHo H T e~ itHa . (3.1) 

Then the solution is: 

W(t) = l-i [ dsH T (s)W(s) 

Jto 

= ! + £)(-*)* / dg i f 1 ds 2- f " 1 ds k H T ( Sl )H T (s 2 )...H T (s k ). (3.2) 

Using the ciclicity of the trace and the definition of W(t) one rewrites Eq. fl2.14p as follows: 

(Jo(t))rrf= lim Tr r{pi A) W*(t)J a (t)W(t)}, J a (t) := e im ° J a e~ UH °. (3.3) 

A— >oo 
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It is clear that by replacing W(t) as given by Eq. (l3.2[) in Eq. (|3.3p one obtains a full expansion of 
the current w.r.t the tunneling Hamiltonian Ht- Our strategy is to show that one can perform the 
thermodynamic limit on each term in this expansion. Let us make a few remarks on the structure 
of these terms and give the main steps we follow for calculating them. 

i) Given the structure of W(t) and Ht the current will be a series of monomials containing 
combinations of creation/annihilation operators from both the leads and the sample. However, 
due to the particular tensor product form of p^ , the particle number conservation requires that 
in all monomials with a non- vanishing contribution to the trace, the number of creation operators 
should equal the number of annihilation operators separately for the sample, and for each lead. 
It also means that each such monomial contains an odd number of Ht's and is of even order in r 
since the current operator itself is proportional with r. 

ii) In order to simlify notation, we write a&(x) instead of a&(\x)), that is we identify the site 
x with the basis vector \x). We deal with the operators acting on jF$ by systematically inserting 
the projections of many-body states {|^)(^|} between any two Ht's- Using the matrix elements 
A vv i introduced above (see Eq. (|2.16p ) and the shorthand notation df(x) — e ltH ° a# (x)e~ ltH ° one 
has for example: 

(v Q ,H T {s k )v k ) =J2e iSk{E "°- EUh) [Av Q v k (m« k K k (O ak ) + <„ fc < (0« J] . (3.4) 

oik 

Note that A vv i couples many-body states whose particle number differ by at most one. Also, A vv i 
does not depend on A, thus the thermodynamic limit is only relevant for terms of the type: 



Tr^ {p[ A) a* 1 (0 Q1 ) . . .a*™ (0 a2N ) } . 



iii) Next we change the representation of the operators using the eigenstates <fq ak 01 the leads' 
Hamiltonian: 

«f fc fc (0 Q J = e^#*^V#* (0 a >#* , (3.5) 
i"k 

where we introduced the notations: 



J# k 



- fOT <*, " V(°aJ for a <W, 



and ip(O a ) = (tpq a ,O a }- The general term on which one should perform the thermodynamic limit 
reads as follows: 

T,T,T,^* iS1£{qai)+ --- +l6 * 2NS2N£{ ^ ( 3 - 7 ) 

S q S g 

where in the trace above there are precisely N creation and N annihilation operators from the leads. 
We introduced the shorthand notations a := (ai, .., «2Ar), qs ■— <Zai > ••, Q_a 2N ancl # := #1 , ■■H Z 2N ■ 
The trace is further calculated using the Wick theorem (see [9]) which holds because the leads are 
noninteracting. The idea behind the Wick procedure is to systematically use the anticommutation 
relations in order to reduce the monomial of order 2N to a sum of monomials of order 2N — 2. 
The simplest case corresponds to all six combinations for N = 2. For example: 

TT ^Ap { L )a qai a* qa2 a qa3 a qai } = -5 qaiqa3 S qa2qa J ai (e qai )f a2 (e qa2 ) 

+ &q ai q a Jq a2 q a3 fa 1 { £ q ai )fa 2 { £ q a2 )i ( 3 -8) 

where we used the cyclicity of the trace, the identity a* Qi — e l3£qa ^ pLa qai and the well known 

fact ^f l {Pl o-q a Q>qp} — S a pf a (e qa ), where f a is the Fermi function associated to lead a. One 
can easily show that all allowed combinations of 4 operators can be expressed in terms of products 
//: // an(1 / /) where / = ! — /. Also, it is important to observe that due to the Kronecker 
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symbols the sums over g's are reduced and one actually obtains products of terms which are of 
the following type: 

J2 e±l[S ~ S ' )eq "M^)(^,^)(^,0p) = (0 /3 ,e ± ^- s '^ A) / /J (/ 1 ( A )) 0/J ), (3.9) 

9/3 



The second term appears from combinations containing / /. 

For terms of higher order one proceeds in a similar way using the general formula (see Eq. (24.36)) 
in Ref.0: 

Tr^ { P { L A) a*\ ..a#f N } = {a*^ , a*^ (e Qai )Tr^ L {p^af^ a*^ ..a*™ } 

- i a l\ > a l\ ( £ ? Q1 ) Tr ^ {p { L )a t 2 2 a t\ -4z } + ■■■ 

Thus we have shown that the thermodynamic limit is to be performed only on factors like in 
(|3.9j) . We give this result as a general lemma: 



Lemma 3.1. LetJ\f A be the set {0, 1, . . . , A} with A < oo. Let hoo be the discrete Laplace operator 
on the halfline Afoo with Dirichlet boundary condition at —1, and h A is the restriction of h^ on 
A/a with Dirichlet conditions at — 1 and A + 1. Let F be any continuous function defined on the 
interval [—2tL,2ti]. Then we have: 



1 f 2tL / E 2 

lim (0, F(h A )0) = (0, Fih^O) = — / 4 / 1 - ^. (3.11) 

a->oo nt L j_ 2tL y 4^ 

Proof. Fix e > 0. The spectrum of all /ia's is contained in [— 2it, 2ti]. Because F is continuous on 
this interval, it can be uniformly approximated with polynomials. The Weierstrass approximation 
theorem says that there exists a polynomial P e (x) — J2j=o a j x ^ such that 

\\F - PeWoo := sup \F(x) - P e (x)\ < e/3. (3.12) 

x€[-2t L ,2t L ] 

The spectral theorem implies that ||-F(^4) — P e (^4)|| = \ \F — P € ||oo for any self-adjoint operator A 
whose spectrum lies in [— 2ij,, 2ti]. Thus we can write: 

|(0,f(U0)- <0, J P e (/i oo )0)j < e/3, \{0,F(h A )0)-{0,P e (h A )0)\<e/3, VA > 1. (3.13) 

It is very important to note that the above estimate holds true uniformly in A. Now let us remark 
that there exists A e sufficiently large such that 

(0,P 6 (/ioo)0) = (0,P e (%)0), VA > A £ . (3.14) 

The explanation is that h A \0) = h^O) if k < A, because we cannot reach the 'other' boundary 
after less than A steps. Thus choosing A c larger than the degree of P e is sufficient to conclude 
that P e (h A )\0) = P e (hoo)\0). Now using flSjJJl and ([3~T4|) we have: 

\{0,F(h A )0) - (0,^(^)0)1 < 2e/3 < e, VA > A e , 

and the proof is over. □ 



After applying the Lebesgue dominated convergence theorem on the iterated integrals of (|3 . 3[) . 
we arrive after some work at a rough estimate of the form: 

j-2n 

\I a AT)\<^J 2n C^— r (3.15) 
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where C is some positive constant. Thus I a ,t{') is entire in t. But this estimate only says that 
the transient current cannot grow faster than an exponential of the type T 2 e Crt , which is not very 
useful if t is large. But at least if r\ is chosen such that t/tj is small enough, then (|1.3|) holds true. 

4 Off-resonant transport 

Before starting our calculations we review the Landauer formula for non-interacting electrons [5] 
which was proved to give the steady-state current both for discrete and continuous models at 
arbitrary bias [TJ 1151 [5] . The reason to make some connection between the Landauer formula and 
our results is twofold. On one hand any calculation in the interacting case should lead to this 
formula when the interaction strength U is set back to zero. On the other hand one can get some 
general facts about the expansion of current in powers of the lead-dot tunneling r. The Landauer 
formula gives the steady-state current in the lead a: 

I«,oo(t) = E n dE (U E ) - UE))\T ai {E)\\ (4.1) 

7 J-2t L 

where the transmittance T ai {E) is defined as follows (see [B]): 

T^(E) = ^Jl-^(m a ,(h s -E-^(f(E)n T ) my), (4.2) 



Trti V 4« \ ' V t L 



where we introduced the orthogonal projection on the contact sites mp, LLp := \m^)(mp\ and 
£ ( z ) = C+ 0) if Im(z) > 0, Ci(z) = C- (z) if Im(z) < 0, where: 

(±(z) = ^ r Ti^l-z 2 /(K)> z ? ((-oo, -2ti]U[2* L ,oo)). (4.3) 

If all eigenvalues e\ of h$ are far away from the spectrum of the leads, then the ergodic current 
becomes analytic near r = and the leading term is of order r 4 and coincides with (|2.22l) . 

In contrast, if some eigenvalue of the sample is inside (— 2^,2^), the ergodic current 
has a completely different behavior with t. For simplicity, assume that all other eigenvalues of 
hs are outside [— 2tL, 2£l]j while e\ is non-degenerate and corresponds to an eigenvector </>, i.e. 
hs<fi = &\4>- Then following [6] one can prove: 



(r) 

r 4 r 2tL 1 - T?r 

E / dE\4>{ mi )m ma )\ 2 (UE)-f,(E))- ^ ,M2 +0 ^> 

„J-2tr. \E- e\+ ttCi (E)(4>, U T 0)r 

(4.4) 



^ in/ + \S^\t ( \ fi M l^( m 7)l 2 l^( m «)l 2 , nn 
r \C(e x ,t L )^[f a (e x )-Me x )] \^ mp )\2 



where C{e x ,ti,) is some constant. It is clear that this expression has nothing in common with 
(|2.15j) . which only contains f a and not differences of Fermi functions. 

4.1 Proof of (ii): Sequential tunneling contribution 

In this Section we calculate the first two contributions to the steady-state current, that is the 
terms of order two and four in the transfer Hamiltonian. Using the identity: 

g-rfffgrf/fo = l-i [ dse~ isH H T e isHo (4.5) 
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and expanding the unitary evolution e ttH up to the 2nd order in Ht one gets from (|3.3p and 

POO nt 

C a , 2 (v)=iV dte^ 1 ds{[H T {-s),J a ]) Te i. (4.6) 
Jo Jo 

By replacing J a and Ht one arrives after straightforward calculation at the following expression: 
[H T (-s), J a }} = r 2 ^^|^ Q (0 Q )| 2 |^ o | 2 (l-/ Q ( £? J( e -(^-^o+^) +C . C ) 

= r 2 ^ ^ |2( e »(^-^o) ( 0a , (1 - f(h L ))O a ) + c.c) 

- r 2 ^|^ o | 2 ( e -^-^o)( 0ct , e -^/(/, L )0 Q )+ C . C ). (4.7) 

In the thermodynamic limit one has: 

{0 a ,e- ishL f(h L )0 a )= I L dE\^(0 a )\ 2 e-' lsE f a (E), (4.8) 

J-2t L 

where ip^ denotes the generalized eigenfunction of the semiinfinite lead corresponding to energy 
E (see (|2.6[1 ). By performing the time integral over s the contribution of order r 2 to the transient 
current is obtained as: 

2r 2 £ dE\^ a )\ 2 (Wol 2 (l ^ f S"^ ~ UlAtJ ^r^ ) , (4.9) 

where for simplicity we omitted to write the energy dependence of the Fermi functions and we 
introduced the notations: 

At (E):=E 1/ -E Uo ±E (4.10) 
Then we perform the final time integral and use the identity: 

limn/ dte-v* = 2tt5(A± J. (4.11) 

to arrive at Eq. (|2.15[) . 

4.2 Proof of (iii): Cotunneling 

The 4th order contribution to the current follows from the expansion of the unitary evolution up 
to the 3rd order in the tunneling operator Ht'- 

T 4 C aA {f])^i f ds x f 1 ds 2 f ds[(e- lSlHo HTe lS2Ho H T e- lS2Ho e* SlHo J a e-< H( >H T e ls 'i Ho ) Icf 
Jo Jo Jo 

-i f dsi f f 2 d S3 (e~ lSlHo HTe lS3Ho HTe- lS3Ho e lS2Ho H T e^ S2Ho e lSlHo J a )^ + c.c. 

Jo Jo Jo 

(4.12) 

In order to achieve a more explicit form of C a ^(rj) we follow the same steps as in the proof of 
the thermodynamic limit, that is we insert the many-body states of H$ in order to deal with the 
operators acting on J\s, then we switch to the proper basis of Hl and finally use the Wick theorem 
for all non- vanishing combinations of the type Tr{a^a^a*a^}. The calculations are tedious but 
straightforward. We find that there are 48 terms contributing to the cotunneling current. At the 
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next step we perform the time integrals. It is sufficient to calculate the real part of this integrals 
because 24 terms are the complex conjugates of the remaining ones. Moreover, one notes that 
there are only two types of integrals: 



nt nsi pt 

B\{t) = I ds\ I ds% I dss cos(six + S2y + S3Z) 
Jo Jo Jo 



(sintz — smt(x + y + z) + smt(x + y)) 



zy(x + y) 



+ (smt(x + z) — sintz + s'mtx), (4-13) 

xyz 

pt rst rs 2 

B2(t) = / dsi / ds 2 / ds3 cos(six + s 2 y + s 3 z) 
Jo Jo Jo 

s'mt(x + y) smtx sinte sint(x + y + z) 



zy(x + y) zx(z + y) xyz z(y + z)(x + y + z) ' 



(4.14) 



where x, y, z contain two many-body energies of H$ and energy of one or two electrons from the 
leads (an example is x — E u » — E UQ + e qi , y — E u > — E u >> — e Q2 , z = E v — E u > + e q2 ). Then one 
has to perform the thermodynamic limit, to calculate the integral over time and take the limit 
?7 — >• 0. This final step brings in plenty of delta functions. Our first off-resonant condition was 
that E v — E v i — e q 7^ if the number of electrons in the MBS \ v') differ by one. Our second 
off- resonant condition implies that E v — E v i ± (e qi + e q2 ) 7^ , for any pair of many-body energies 
E V ,E V whose particle numbers differ by two. By analyzing all combinations of x,y,z it follows 
that the remaining off-resonant terms arise from 5{x)/ zy for B\ and from S(x + y)/zy for B 2 . In 
these terms the delta functions impose conditions of the form E v * — E VQ + e qi — e q2 = 0, which 
means that the dot initially in the state \vq) passes to the state \v') by exchanging two electrons 
with the leads. This process is called cotunneling in the physical literature, because the electrons 
now tunnel pairwise. After collecting all these terms and taking advantage of some cancelations 
one arrives at the final expression for the cotunneling current given by Eq. (|2.17|) of the theorem. 

Let us make a few remarks on the cotunneling current. From the sequence of A's appearing 
Eq. (|2.18|) one observes that the cotunneling processes always imply different leads. Take for 
example the 3rd term. It describes the following sequence: an electron with energy E enters the 
dot from the lead a, while the second electron of energy E' = E — E v * + E Uo leaves the dot to lead 
7. The remaining two terms described the reverse process: the electron tunnels back from the lead 
7 and the second one tunnels out to lead a. The other terms can be described in a similar way. 
Also note that the cotunneling contributions cxplicitely contain the initial state of sample vq. 

A natural question is what we can say about the cotunneling current in the non-interacting 
case. Let us recall here that e\ are the eigenvalues of /is, i.e hs4>\ — e\4>\. Then the operators 

(|ro a )) and a^(|m 7 )) appearing in the coefficients A in Eq. (|2.17[) can be written in terms of 

:— a^(\(f>\)). Moreover, the sums over the many-body states of H$ allow one to recover the 
resolvent (Hs — E VQ — e)^ 1 and also the Fermi-Dirac operator / Q!7 (i?s — E VQ — e). As an example 
we consider the 2nd term in Eq. (|2.18l Introducing the notation f-y(E) := xl(E)J 1 (E) one has: 

M 2 :=- £ ^Jd^M^. A* V0V (m a )A vv . (m^A*^,, (m 7 )A^ (m a ) 

= ^2 (^^^ m a)(m.y,(j)x 2 )((f>x 3 ,m 1 )(m a ,(l)x i )(l- f a (E))f 7 (E + E^ - E Vo ) 

Ai ,A2,A3 ,A4 

x (vo,a* Xl (Hs - E V0 - E^axJ^Hs - E VQ - E)a* X3 (H s - E VQ - E)- x )a M v Q ) 

ST I A \l A MA V A \ ( X ~ fa)M E + e >2 - eAi) 

= 2^ { < P\i, m a){m 1 ,(px 2 ){(j)x 3 ,rn 1 ){m a ,(f>xJ 7^73 )[e~^ ) 

A\\\ \ A4 / \ Ai / 

1 ,A2 ,A3 ,A4 

xrT ^s{Psa*x 1 ax 2 a*x 3 a Xi }- (4.15) 
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In the above calculations we used pull-through identities like (Hs — z)~ 1 a\ = a\(Hs — z — e\)~ x 
or a* x f(Ho — z) = f(Ho — z — e\)a* x . Now the only thing we should do is to use the Wick theorem 
for the trace in the last line (the theorem now holds as the interaction is absent): 

Tr^ s {psa^aA 2 a^ 3 aA 4 } = n Xl (6x^X3X^X3 + S Xl \J\ 2 \ 3 (1 - n\ 2 )), (4.16) 

where n x = Tr^{psa* x ax} and 1 — n x = Trjr{p s a\al}. 

The remaining terms in Eq. d2.18p have to be manipulated in the same manner. Collecting all 
of them one observes that all products of Fermi functions vanish, and all factors like / will only 
appear as f(E) with E £ [—2^,2^] which allows us to drop \l- Then the cotunneling current 
takes the following form (which does not depend on the particle number n x ): 

r 4 r 2tL 

r 4 C a>4 (0 + ) = ^- 3 -X] / [l - E 2 /(4t 2 L )]\(m a ,(hs - Ey'm^iUiE) - f^E)) dE. (4.17) 

n l L 7 J-2t L 

One recognizes at once the 1st term in the expansion of the Landauer formula (14. ip w.r.t t in 
the off-resonant case. So as expected, the off-resonant transport is still described by a Landauer 
formula in the non-interacting case. As expected, in this case the steady-state current does not 
depend on the initial state of the sample. 

5 Numerical simulations of the transient regime 

Let us consider the same two-site system as the one in Remark l2.5l As we have already mentioned, 
one can numerically compute transients via the generalised Master equation (GME) method [13]. 
The main idea behind this method is to write down an equation for the reduced density operator 
(RDO) p r (t) := Tijr L {p}. Note that p r (t) only acts in the Fock space of the sample. Its derivative 
w.r.t time gives the evolution of the particle number in the sample which in turn is related to 
the currents flowing to and from the leads via the continuity equation. The method is usually 
formulated in terms of Liouvillians (see e.g. |18) for relevant equations). Although the main regime 
considered in other papers is the resonant one, here we pay more attention to the off-resonant 
regime. We are motivated by the fact that in our paper the transient current due to sequential 
tunneling processes is given by a rather simple analytical formula Eq. (|4.9l) . which should be the 
main contribution on a time scale of order 1/r. 

Moreover, since GME also works for the resonant case, it would be a proper tool to compare 
the two-regimes. The off-resonant setup is achieved by taking a small hopping constant on the 
leads and by globally shifting the leads' spectrum ctQil) = [— 1th + -Eshift, 2£l + £7 s hift]- The bias 
window [ptR, pl] is also fixed such that all the many-body states of the sample are below it. The 
time-dependent currents in the left (L) and right (R) leads are presented in Fig. la. 

The convention for the sign of the currents is as follows: Jl is positive if it flows from the left 
lead towards the sample and Jr > if the current flows from the sample to the right lead. The 
steady-state regime thus implies Jz,(i) = Jr(^) for some t. Instead of this one notices that both 
currents exhibit modulated oscillations around zero and no steady-state is achieved, although the 
amplitude of the oscillations decreases in time. This behavior could be predicted by our analytical 
result (see Eq. (|4.9[) ). However, if one performs the ergodic limit the results converges to zero in 
the long-time limit, as clearly seen in Fig. lb. 

The transport in the resonant regime is shown in Fig. 2a for two initial conditions of the isolated 
quantum dot \i/q) = 1 10) and \vo) — 1 00) . In this case we consider a larger th and the bias window 
is chosen such that the first state of the dot is below it while the other ones within the bias window. 
Notice that in this case the parameters are set such that <j{h£) covers the entire spectrum of hs- 
The transients are quite smooth and the steady state is achieved around t ~ 225. In this case 
there is no need to consider the ergodic limit. 

Remark 5.1. In the resonant regime, the steady-state current does not depend on the initial 
condition of the sample. This has already been rigorously established both in the non-interacting 
case [U [8] and for weakly interacting systems [7] . 
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Figure 1: (Color online) (a) The total transient currents Jl and Jr as a function of time in the 
off-resonant regime, (b) The 'ergodic' currents. Other parameters: U — 0.5, r = 0.5, II = 0.1, 
E s hih = 6, i±l = 7, = 6. 
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Figure 2: (Color online) (a) The total transient currents Jl and Jr as a function of time in the 
resonant regime. Two initial conditions were considered N = 1 corresponding to one electron 
on the lowest state and N = corresponding to an empty sample, (b) The occupation of the 
many-body states with n electrons and the total occupation. Other parameters: U — 0.5, r = 0.5, 
th = 1.5, -Eshift = 3, fiL = 5, p,R = 2. 



Remark 5.2. If the quantum dot is initially empty, the current on the right lead starts by being 
negative, which means that this lead actually feeds as well the dot. Fig. 2b shows the charge 
that accumulates in time on the many-body states containing n particles, and the total charge 
n tot (the curves correspond to the initial condition \uq) = |10)). The reading of the numerical 
results is straightforward. The single particle state are depleted in favor of the two-particle state 
|11). In the steady-state regime the latter contains in average one electron, because the state |11) 
contained within the bias window charges/discharges by back-and- forth tunneling of one electron 
from the leads. 

Remark 5.3. The results presented in this section were obtained by numerically implementing 
and solving the integro-differential equation for the reduced density operator which served us to 
calculate the transients. A legitimate question is how one could use the GME method if interested 
only in the steady-state regime? The most tempting step is to assume that a steady-state exists, 
which in terms of the RDO means that lim t _ i . 00 p r (t) = 0. If so, then one can calculate the 
stationary RDO from the GME equation and derive the steady-state currents. This strategy is 
extensively used in the physical literature. Our analysis shows that in the off-resonant regime such 
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an approach is not justified because there is no steady-state. The correct procedure is to work out 
the time-dependent equations and calculate various contributions to the ergodic current which is 
the meaningful quantity to look at. 

6 Conclusions 

We have presented a rigorous approach to the cotunneling transport in weakly-coupled interacting 
quantum dots. Using the expansion of the transient current in powers of the lead-dot coupling 
parameter r we analysed the leading order contribution (i.e. C(r 4 )) of the ergodic current which 
is the relevant quantity to consider in this regime. Explicit calculations for elastic and inelastic 
cotunneling contributions to transport were presented. For non-interacting electrons one recovers 
the Landauer formula. For a simple two-level system, we show that in the interacting case the 
cotunneling current depends on the initial many-body configuration of the dot. To our best 
knowledge, this memory-effect has not been reported before. An explicit formula for the ergodic 
sequential tunneling current (i.e. C(r 2 )) is given. This contribution vanishes in the cotunneling 
regime but the transient sequential tunneling does not reach a stationary state. These results 
are also recovered through numerical simulations via the generalized Master equation method. 
This method allows calculation of transient sequential tunneling currents. A generalized Master 
equation containing higher order terms has been recently reported [19J.[11J. This motivates a 
thorough rigorous analysis on the existence of a stationary regime for the reduced density operator. 
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